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Abstract 

Based on the phosphorelay kinetics operative within BvgAS two component system we propose a mathematical framework for 
signal transduction and gene regulation of phenotypic phases in Bordetella pertussis. The proposed model identifies a novel 
mechanism of transcriptional interference between two promoters present in the bvg locus. To understand the system behavior 
under elevated temperature, the developed model has been studied in two different ways. First, a quasi-steady state analysis has 
been carried out for the two component system, comprising of sensor BvgS and response regulator BvgA. The quasi-steady state 
analysis reveals temperature induced sharp molecular switch, leading to amplification in the output of BvgA. Accumulation of a 
large pool of BvgA thus results into differential regulation of the downstream genes, including the gene encoding toxin. Numerical 
integration of the full network kinetics is then carried out to explore time dependent behavior of different system components, that 
; qualitatively capture the essential features of experimental results performed in vivo. Furthermore, the developed model has been 
utilized to study mutants that are impaired in their ability to phosphorylate the transcription factor, BvgA, of the signaling network. 
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1. Introduction 
I ' 

I ^ I One of the important functional aspects of living organisms 
is to respond to the sudden changes made in their environment, 
and to make appropriate changes in the cellular or subcellu- 
lar level for survival. Direct manifestations of such changes 
at the subcellular level are the expression/repression of single 
or multiple genes controlling different functional behavior of 
an organism (Alon, 2007). To achieve this, living system uti 



As a response to temperature elevation in the environment, 
the response regulator BvgA becomes active (the phosphory- 
lated dimer) within each bacterium, which in turn exerts a posi- 
tive feedback on its own operon, the bvg operon. Positive feed- 
back loop thus increases the active form of BvgA in a switch 
like manner In other words, once the BvgAS two-component 
machinery becomes operative, large pool of active BvgA ei- 
ther repress and/or express several downstream genes where the 
phosphorylated dimer of Bv gA plays the l eading role by acting 



o ' ' / ^ c J phosphorylated dimer or Bv gA plays the l eadmg role by actmj 

^ lizes concerted bio chemical ne twork composed of sever alfegd- transcription factor (TF) (Tsteffen et allH^ . In the labora 



back mechanism ( Tyson et al. , 2003 : Tyson and Novak , 2010l) 
The human pathogen Bordetella pertussis, a gram negative 
bacteria and causativ e agent for the disease whooping cough 
tlPreston et al. , 20041) . is no exception to the aforesaid behav- 
ior. At 25 °C, while freely moving in the environment their 
pathogenic properties remain dormant. But, when they are 
.within the host at 37 °C, their virulent properties come into 
play. In the laboratory the reverse effect, i.e., suppression of 
pathogenic behavior is observed using MgS04 o r nicotinic acid 
(iBeier and GrossL l2008t ICotter and Jones . l2003h . The virulent 
behavior of B. pertussis within host, in response to sudden en- 
vironmental change, has been experimentally studied and has 
been found to be operative through BvgAS two componen t 
system (TCS) dBeier and Gross! bOOsUCotter and Jonesi (20031) . 
The TCS comprises of transmembrane sensor BvgS and re- 
sponse regulator BvgA where signal flows through this pair via 
a four step (His-Asp-His-Asp) phosphorelay mechanism. 
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tory condition at 37 °C and in the absence of MgS04 or nico- 
tinic acid BvgA activates transcription of virulence activated 
genes {vag), as well as represses transcription of virulence re- 
pressed genes iyrg). Due to this reason bvg locus was earher 
called vir due to its connection to virulence dBeier and Gross! 
2008t ICotter and Jonesl Eool IWeiss and Falkowl Il984 . In B. 
pertussis, vrg loci encodes outer membrane whereas in B. bron- 
chiseptica, vrg controlled genes encode motility and survival 
from nutrient limitation condition. In Bordetella spp., vag loci 
encodes genes responsible for adherence, toxins (including per- 
tussis toxin in B. pertussis, a type III secretion system and Bv- 
gAS itself (due to autoregulation). 

Based on the binding affinity of TF to the respective promot- 
ers, different types of downstream genes are regulated in Bor- 
detella spp. and have been broadly grouped into four classes, 
e.g., class 1, class 2, class 3 and class 4 dBeier and Gross . 
2008tlCotter and Jonesll2003h . Class 1 genes encompass genes 
that are responsible for encoding toxins, such as adenylate cy- 
clase {cyaA-E) and pertussis toxin iptxA-E). Class 2 genes ex- 
press proteins responsible for adherence, such as fliaB that en- 
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codes filamentous hemagglutinin. Among all the four classes of 
genes, class 3 genes show a unique beh avior, although its func- 



tional activity is not kn own till date jBeier and Grossl 12008 



Cotter and Jonesl 120031) . The only well characterized class 3 



gene found in B. pertussis is known as bipA. The final one, class 
4 genes have been reported to encode frlAB in B. bronchisep- 
tica and is responsible for motility. It is important to mention 
that expressions of class 3 and class 4 gene are not observed in 
B. pertussis under the influence of temperature elevation. To be 
specific, class 3 gene expression has been observed in B. per- 
tussis only under the influence of intermediate concentration of 
MgS04 and class 4 gene under low concentration of MgS04 
in B. bronchiseptica jBeier and Grossl l2008t ICotter and Jones , 
20031). In terms of vir regulated genes class 4 gene thus be- 



longs to vrg whereas class 1 and class 2 genes belong to vag. 
Expression and/or repression of the four classes of downstream 
genes is controlled by strong and/or weak binding sites (for TF) 
present in the promoter region of the respective genes. Among 
these, promoter region of class 4 gene has the strongest aflin- 
ity for TF. Promoter region of class 2 and class 3 genes have 
medium affinity for TF, whereas promoter region of class 1 gene 
has the weakest affinity for TF. On the basis of the promoter re- 
gions' affinity for TF it is thus expected that expression and/or 
repression of four classes of downstream genes in Bordetella 
spp. would show a differential pattern in their temporal dynam- 
ics. 

Keeping these aforesaid phenomenological information in 
mind we have developed a mathematical model based on bio- 
chemical interactions taking place within B. pertussis under the 
influence of temperature elevation. The objective of present 
work is twofold. First, we aim to understand the molecular 
switch operative in BvgAS TCS and to identify the key players 
responsible for amplification of TF. Second, through our model 
we aim to regenerate qualitative features of the network and to 
mimic different phenotypic states of B. pertussis under temper- 
ature elevation, as well as their expression level due to different 
mutation. 

2. The Model 

To understand the mechanism for temperature induced ac- 
tivation of bvg locus and differential regulation of the down- 
stream genes, we propose a kinetic model in the following. 

2.1. The bvg locus 

Experimental studies in B. pertussis suggest multi-promoter 
activities in bvg operon (Ro v et al. . 119901: ] IScarlato et al.l . Il990l 
). Out of the four promoters, Pas\, Pasi, Pas3 and Pas4, 
present in the bvg locus (see Figure [U, only Pasi is known 
to be constitutively active under non-inducing condition (25 
°C) and is bvg independent. After induction (37 °C), activity 
of the Pas 2 promoter goes down while the other three pro- 
moters (Pasi, Pas 3 and ^"^54) become active. As shown in 
Scarlato et al.l (Il99ll) . at 37 °C, Pas 1 shows maximal level of 
activity compared to ^"^53 and is on within < 10 minutes of 
induction. The amount of transcripts generated from Pas 3 is 
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Figure 1 : (color online) Schematic presentation of hvg locus and signal trans- 
duction in BvgAS two component system. The dashed Hne presents the feed- 
back by phosphorylated dimer of BvgA on its own operon. The dotted line is 
for the production of dimers of BvgS and BvgA. For simplicity the mRNAs are 
not shown in the diagram. 



very low and have been reported to be hardly detectable. The 
Pas4 promoter shows same level of activity as Pasi but pro- 
duces antisense RNA. Although activity of Pas4 promoter and 
its product, the antisense-RNA, is known, the target of the an- 
tisense RNA is not known till date. In passing it is important 
to mention that multi -promoter activity in the operon of TCS as 
observed in B. pertussis, has a lso been obse rved in other human 



pathogens (iChauhan and Tvag i, 2008; Do na et all 120081) . 



To study functioning of the bvg locus we consider only the 
activity of two the promoters Pasi and Pasi in our model, as 
reasonable amount of experim ental data is available in the liter- 
ature for these two promoters ( Scarlato et al. 199ll) . Both these 
promoters are typical example of tandem promoter, containing 
a conserved region of a; 10 base pairs between upstream of 
Pasi and transcriptional start site (TSP) of Pasi- In the model, 
we designate the constitutive form of Pas 2 under non-inducing 
condition as Pasic- Once induced, TF interacts with both Pas\ 
and Pas2 and makes them active 



PaS2c+A2P ^ PAS2a, 



PAS\i+A2P 



AS\a- 



(1) 

(2) 



In the above equations PAS2a is the active form of Pas2\ and 
Pasi I and Pasxu are inactive and active form of Pasi pro- 
moter, respectively. Due to the presence of conserved region of 
^10 base pairs, RNA polymerase (RNAP) for Pasi traversing 
through downstream of Pasi now interferes with the binding 
of TF (and RNAP for Pasi) to upstre am of Pas2 causing tran- 
scriptional interferenc e (see Figure |2) ( Buetti-Dinh et al. , 20091 : 
Shearwin et al.l l2005h . During this process Pas2 and Pasi act 
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Figure 2: (color online) Schematic presentation of transcriptional interference 
between Pas I and PaS2- RNA polymerase (red blob) starts its journey from up- 
stream of Pas 1 promoter executing initiation, elongation and termination (fol- 
low the dotted arrowhead). During termination it interferes with the TF (A2p, 
magenta dimers with yellow blob on top), causing transcriptional interference 
and downregulation of Pas 2 promoter activity. The vertical solid aiTowheads 
presents binding/unbinding process between TF and Pas2 promoter site. The 
filled bar at the bottom is for overlapping » 10 base pair region between up- 
stream of Pas I and TSP of Pas2- 



as sensitive and aggressive promoter, respectively. Although a 
detailed kinetic mechanism of transcription al interference has 



been proposed and verified experimentally (iBuetti-Dinh et al 



2009h . we use the following notion to keep the model simple 



Pas 2a ^ Pasii, 



(3) 



where the information of upstream inhibition are put together 
in the rate constants. The rate constant ka in Eq. ([Sj contains 
the information of RNAP coming from Pas i causing transcrip- 
tional interference at Pas2, 

PAS2a+ RNAP %P ASH. 

Considering the pool of RNAP to be constant one can absorb 
it into the overall rate of the interference mechanism and write 
RNAP k'.2 - ka, the overall rate constant of the interference 
process given in Eq. After causing the interference the 
RNAP coming from Pasi falls off from Pas2 thus giving chance 
to the later to back to the active form again, which is modeled 
as the backward reaction with rate constant ka2 in Eq. (O. This 
helps the Pas 2 promoter to maintain a low activity even after 2 
hours of ind uction as evide nt from the experimental data (see 
Figure 3B of IScarlato et al.l (11991.) and Figure |6]of the present 
work). To keep track of the transc ripts g enerated due to Pasi 



and Pas2, following lScarlato et al.l(ll99 1). we consider two iso 



forms of mRNA generated from the bvg locus, niAsi and mAS2, 
respectively. At 25 °C, mAS2 is constitutively produced from 
the constitutive state of Pas2 promoter 



AS2c 



mAS2- 



(4a) 



After 2 hours of induction, level of mAS2 goes down but 
still maintain s a low level of expression (see Figure 3B of 
Scarlato et al. d 1991 1)-) which we assign to the residual activity 
of PAS2a- After 20 minutes of induction, production of mAS2 
remains still on, reaches a maxima, and then goes down which 
we attribute to the suppression of PAS2a and the formation of 
PAS2i altogether. Thus after induction. 



Pas 2a — * mAs2, 



Unlike the transcripts generated due to Pas 2 activity, generation 
of niAs 1 is solely governed by the active form of Pas i promoter 



Pas la — > niAs i ■ 



(5) 



Finally, we consider natural degradation of both the transcripts 
generated from Pas\ and Pas2, 



kd.,„ kjj„ 

niASi — > 0,mAS2 — > 0- 



(6) 



2.2. The two component system 

Once transcribed from the locus, we consider translation of 
two isoforms niAsi and mAS2 into dimers of sensor, BvgS (5 2), 
and response regulator, BvgA (A2) proteins. 



niAs 1 — > niAs 1+^2 

k,. 

mAS2 — 



mAsi 



mAS2 



mAS2 +S2, 
mAsi +A2, 
mAS2 +A2. 



(7a) 
(7b) 
(7c) 
(7d) 



In reality the sensor and response regulator proteins are first 
translated as monomers and then dimerize (Beier and GrossI 
,2008.) . Since dimers of BvgA, not the monomers, act as tran- 
scription factors for the activation of bvg locus, we have omitted 
kinetics of monomer formation and subsequent dimerization in 
our model and work instead with 5 2 and A2. As we show in the 
following, this simplification does not affect our analysis and 
modeling of the signal transduction network. 

In bacterial TCS, autophosphorylation occurs at histidine 
residue of the sensor kinase, that serves as source of phosphate 
group and transfers the same to the asparta t e residue of re s ponse 



re gulator, acting as sink (Applebv et al.L 11996 : iHochi 



Laub and Goulian , 2007; Mitrophanov and Groismani 



2000; 



2008 



Sto ck et al.L l2000h . This orthodox two-step His-Asp phos- 



photransfer becomes complicated in B. pertussis where signal 
transduction takes place via a unortho dox four-step His-Asp - 



His-Asp phosphotransfer mechanism (lUhl and MiUerl Il996l) . 
where the first three steps (His-Asp-His) take place within 
the sensor protein BvgS and in the last step (His-Asp) phos- 
phate group flows from the transmembrane sensor BvgS to 
the cytoplasmic response regulator BvgA. Mathematical mod- 
eling of two-step phosphorelay has been reported in differ 
ent context of bacterial signal transduction ('Banik et al.', '2009 
Batchelor and Goulian, 2003; Kato et all 12007: Kierzek et al. 



2010; Kr emhnget all 12004 Ishinar et al.n2007t ISureka et al 
2008) . Whereas mathematical modeling of four-step phospho- 



(4b) 



relay has been extensively used for the gram positive bacte- 
ria Bacillus subtilis to understand the sporulation initiation, 
a detaile d account of whic h can be found in the recent re- 
view by iLiebal et al. I dioiol) . In B. subtilis the external signal 
is sensed by a family of five histidine kinase KinA-E and fi- 
nally transferred to the response regulator SpoOA via SpoOF 
and SpoOB. Phospho rylated SpoOA acts a s a TF by control- 
ling over 500 genes dFawcett et all I2OO0I) which are broadly 
classified into two categories by the affinity of SpoOA to their 
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target genes ( Fuiita et"an. 2005|). In pass ing we would like to 
mention the work bv lKim and where a comparative 

study of two-step versus four-step phosphorelay has been un- 
dertaken. In the present study we have adopted a sim plified ap- 
proach of what is proposed bv lKim and Chol(l2006h . Instead of 
detailed three-step (His-Asp-His) phosphotransfer mechanism 
within BvgS we consider autophosphorylation at 5 2, the dimer 
of BvgS, 



^ Sip- 



(8) 



While writing the above equation we have put the information 
of ATP and its interaction with the sensor protein in the rate 
constant kp^siis)'. In addition it is a function of input signal s 
which may be temperature or salt concentration for B. pertus- 
sis. In absence of any signal, i = 0, autophosphorylation at the 
histidine residue becomes nonfunctional, i.e., kj, ,2(0) - 0. At 
this point it is important to mention that exact mechanism for 
the activation of BvgS in B. pertussis under temperature induc- 
tion is not clear from the literature. To incorporate sensing of 
the external stimulus and subsequent activation of the TCS we 
have adopted the mechanism given in |8] Once phosphorylated 
the sensor protein transfers the phosphate group to their cognate 
response regulator A2, the dimer of BvgA, mimicking the last 
step (His-Asp) of the four-step phosphotransfer mechanism 



S2P+A2 ^ S2P-A2 > ^2 + A2P, 

k,j, 



(9) 



where S2P ■ A2 is the Michaelis complex formed by S2P and 
A2. In addition to their kinase activity the sensor protein 
also exhibits phosphatase activity by remo ving the phosphate 



group from their cognate partner (|B atchelor and GoulianI 12003 



Laub and Gouhan , I2OO7I: IShinar et al.l 12007b thus showing bi- 
functional behavior 



kp.f kp^ 
S2+A2P ^ S2-A2P - 

kp.h 



■S2+A2. 



(10) 



Likewise the mRNA transcribed from the bvg locus, we con- 
sider natural degradation of different forms of sensor and re- 
sponse regulator proteins 



^d.p ^d.p ^d.p 

1,S2P — > 0,A2 — > 0,A2P — 



0. 



(11) 



The phosphorylated dimer of the response regulator then exerts 
a positive feedback in the bvg locus thus activating the promot- 
ers of the locus (see Eqs. (1-2)). In addition they control the 
expression and/or repression of several downstream genes. 

2.3. Differential gene regulation 

To model regulation of downstream genes under induced 
condition we use the following TF binding properties of differ- 
ent promoters of downstream genes. The promoters for class 1 
gene contain low aflinity BvgA binding site far upstream of the 
TSP, as a result high level of A2P is necessary to activate these 
genes and they are expressed quite late compared to class 2 and 
class 3 genes. Class 2 promoters contain high affinity BvgA 
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Figure 3: (color online) Phenotypic gene regulation in B. pertussis. The red 
blob is for RNA polymerase and the magenta dimers with yellow blob on top 
are for phosphorylated dimers of BvgA. 



binding site close to TSP and fairly low level of A2P is suffi- 
cient to activate these genes. As a result, expression of class 2 
genes becomes visible within very short period after induction. 
Class 3 gene, bipA, contains high affinity as well as low affin- 
ity BvgA binding site just upstream and downstream of TSP, 
respectively. Once induced, bipA becomes active almost at the 
same time like ftiaB (class 2 gene) with the help of low amount 
of A2P. At a critical concentration of BvgA gene expression be- 
comes maximum and then it starts falling due to BvgA binding 
to the downstream low affinity binding site. The frlAB promoter 
(class 4) has been found to contai n distinguishable Bv gA bind- 



ing site that overlaps with TSP ( Akerlev et al.l 119951) . Thus a 



very low level of BvgA may be able to suppress class 4 genes. 

In Figure|3]we schematically show regulation of four classes 
of genes as BvgA level increases. The binding kinetics of phos- 
phorylated dimer of BvgA {A2p) to the promoters of these four 
classes of genes are given in the Supplementary data. While ac- 
tive, four classes of genes starts transcribing their correspond- 
ing mRNA 



(12) 



where j = 1-4. Likewise the transcripts generated from the 
bvg locus, we consider natural degradation of the four different 
classes of transcripts 



kd.,„ 

mcij — > 0- 



(13) 



Before proceeding further we would like to mention that all the 
relevant symbols designating biochemical species used in this 
section are Usted in Table 1 . 



3. Results and Discussions 

To check the validity of our proposed model, the kinetic 
network developed in the previous section has been translated 
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into sets of coupled nonlinear ordinary differential equations 
(ODEs). To understand functionality of temperature induced 
switch operative in the proposed model we first analyze the ki- 
netics for bv ^ locus and the TCS us ing quasi-steady state ap- 
proximation ( Batchelor and Gouliani 12003). The full network 
kinetics is then numerically integrated for large set of parame- 
ters of the pro posed model and com pared with in vivo experi- 



mental results (IScarlato et al.Lll991l) . 



3.1. Steady state analysis 

According to the model described in the previous section to- 
tal amount of sensor and response regulator proteins can be ex- 
pressed by the following relations 

[St] = [S2] + [S2P] + [S2P-A2] + [S2-A2PI (14) 
[At] = [A2] + [A2P] + [S2P-A2] + [S2-A2P], (15) 

where [A7-]/[57-] ^ ksa,i/kss,i ~ 6 (see Table 2). Thus one 
can express essential features of the dynamics in terms of the 
response regulator protein. To realize the nature of amplifi- 
cation in gene expression at 37 °C and to analyze the steady 
state dynamics we divide the TCS signaling network into two 
modules, the autoregulation module and the phosphorylation 
(or post-translational) module. 

Before proceeding further we would like to comment on the 
input-output relation in the present model. For this, we define 
the influx and outflux of phosphate group in the network as 

Ji ^ kp^s2is)[S 2] and Jo = kp^a2[S2 ■ A2p], 



respectively (IShinar et al.l 120071) . Using the steady state ex 



pression of [5 2 • A2p] (see Eq. (35) of the Appendix 5.4) and 
considering 7,- - Jo at steady state, one arrives at [A^p] - 
kp,s2{s)KMp I kp^a2- Thus, the steady state output of TF in the 
present model is dependent on the input signal and the phos- 
phatase activity of the sensor kinase on the response regulator 
For further analysis we have removed the superscript ss from 
the steady state expressions for notational simplicity. 

3.1.1. Autoregulation module 

Autoregulation module takes care of synthesis and degrada- 
tion of TF mediated by positive feedback loop operative within 
bvg operon. This ultimately leads into total amount of response 
regulator, Aj-, expressed in terms of the TF, A2P, in the steady 
state. To derive such expression we take time derivative of 
Eq. (fTSl l and impose quasi-steady state condition on [S2P ■ A2] 
and [5 2 • A2p] (see Appendix 5.4) thus leading into 



^ l3xF2{A2p) + I32F2{A2P)[A2P]2 

+PiPiiA2p)-kd^p[AT], 



(16) 



which finally yields the desired functional relationship between 

At and A2P 



[At] « l3xF2{A2p) + I32F2{A2P)[A2P]2 
+/33Fi(A2p), 



(17) 



where y6,- = $i/kjp (i = 1,2,3). The first two terms on the 
right hand side of Eq. (fTTl i arise due to Pas2 promoter ac- 
tivity whereas the third term is solely due to Pasi promoter 
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Figure 4: The steady state behavior of the TCS circuit in terms of normalized 
[A2/>] as a function of normalized [At]. Normalization was done with respect to 
the total amount of response regulator at steady state. The solid curves are from 
phosphorylation module, Eq. 12 IK at low and high temperature. The dotted 
switch like curve is due to autoregulation module, (Eq. U7t ). While going 
from low to high temperature, a small change in the [At] value makes a large 
amplification (the solid line with double headed arrow) in the [A2p] value (see 
the discussion in the main text). The leftmost vertical dotted line is for basal 
expression. 



under inducing condition. In the limit of zero phosphoryla- 
tion (/tp,j2 = 0), i.e., at 25 °C, Eq. ([17]) leads to [At] = ySi 
(= ksa,2ktp,2Q I kd,mkd,p)- Thus, under zero stimulus system dy- 
namics is solely governed by basal transcription (ktp^o) and 
translation (^.,0,2) processes. 

3.1.2. Phosphorylation module 

Considering only the phosphotransfer kinetics (see Eqs. ([8} 
[Tol l) and using the relations given in the Eq. (35) of Ap- 
pendix 5.4 we have at steady state 



kt al 

kp,s2[S2] = kdp.s2[52p] - ■T^[S2P][A2], 
J^Mr 

kn a2 kt o2 

-^[S2][A2P] = -F[52P][A2]. 
l^Mp J^Mt 



(18) 
(19) 



While deriving the above two expressions we have again im- 
posed quasi-steady state condition on the two Michaelis inter- 
mediates [52P ■ A2] and [5 2 ■A2p]. After some algebra Eqs. ([TSl - 
[T9] ) provide 



c, 



A2P A2 



(20) 



where Cp = kp^,2(s)KMp/kp^a2 and C, = kdp,s2KMr/k,,a2- Now, 
for [At-] «i [A2] + [A2p], Eq. ([20b gets transformed into 

[A2p] ^ \{C, + Cp + [At]) 

~.J{Cr + Cp + [AT]f -4Cp[ATl (21) 

which is valid for A2P ^ At- Thus steady state output of A2P 
depends on both A 7 and the source of autophosphorylation. 
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kp,s2(s)- For A7 » C, + Cp, Eq. dlTT) yields A2P ~ Cp. Since 
Cp is a function of kp^siis), for fixed set of other parameters of 
the model, a change in kp^siis) brings in a change in the value 
of A2P. In addition, A2P level does not depend on Aj anymore 
but remains dependent on ATP level through the rate constant 

^/),.s2(^)- 

In Figure |4] we show functional relation between 
no rmalized \A2p\ and normaUzed [A^] at steady state 



dMivashiro and Goulianl l2008h . Normalization was done with 
respect to the total amount of response regulator (phosphory- 
lated and unphosphorylated) at steady state. At steady state, 
depending on the level of input stimulus (here temperature of 
the environment), total amount of BvgA protein, Aj, exists 
either in phosphorylated form (Aap) or in un-phosphorylated 
form (A2) of which only the phosphorylated form acts as 
TF. The level of Aip at steady state is controlled by both the 
external temperature and the total amount of available BvgA 
protein. At low temperature only a small amount of Aj gets 
transformed into A2P which increases as the temperature of 
the environment is raised. Enhancement of A2P level out of 
the total At pool due to temperature elevation is shown by the 
two sigmoidal solid curves in Figure |4] using Eq. (l2ll . It is 
interesting to note that even when a large pool of Aj is available 
to be phosphorylated only a small amount gets transformed 
into A2P at low temperature due to low autophosphorylation 
rate (kp^si) of Bvgs. But, at high temperature high kp^s2 value 
changes the scenario by increasing the A2P level (see the 
discussion on amplification in subsection 3.2). 

Formation of pool of A2P due to temperature elevation does 
not work due to phosphotransfer mechanism only. It works 
hand in hand with the autoregulation module as well. As men- 
tioned earlier, in the absence of external stimulus (kp^si - 0) 
the only form of BvgA available is A2 that leads to Ap - P\. 
This gives the basal expression level (the left most vertical dot- 
ted line) in Figure 4. As the system gets switched on {kp^si + 0) 
A2P gets generated due to phosphotransfer and provides a pos- 
itive feedback on the hvg operon. Positive feedback enhances 
the production of A2 which are ready to be phosphorylated at 
a fixed stimulus. This essentially increases the amount of A2P 
out of Aj. This functional relation between A2P and Aj due to 
autoregulation is given by Eq. ( fTTl i and is plotted in Figure 4 
for nonzero kp^^i (the dotted switch like curve). The intersec- 
tion of the dotted curve and the solid curve at a particular tem- 
perature (low or high) gives the maximum limit of available 
Aj that are available to be phosphorylated. Although, the dif- 
ference between the intersection at low and high temperature 
shows that due to small change in A 7- there is a huge change in 
available A2p. Thus, the intersections of the dotted curve due 
to autoregulation module and the solid lines due to phospho- 
rylation module give a qualitative idea of amplification in the 
steady state output of A2f when temperature of the surrounding 
is increased. 

3.2. Amplification 

The aforesaid discussion gives an idea of how the two mod- 
ules work together to generate the pool of A2P out of the total 
pool of BvgA. Through western blot technique iPrugnola et al. 
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Figure 5: Semilog plot of amplification factor, A, as a function of input stimu- 
lus, kp^sl- 



gave a qualitative picture of amplification of the BvgA 
protein due to temperature sh i ft. Gu ided by the experimental 



observation of IPrugnola et al. ( 19951) and considering the am- 



plification is due to the accumulation of large amount of A2P 
within the cell (see the preceding discussion and Figure (|4)), we 
now look at the response of the system in terms of accumulation 
of A2P as a function of external stimulus. To this end we use the 



concept of amplification factor. A, proposed by Koshlandetal 
( Il982h . In our notation A, can be defined as 



AA2P/Ai 



{^2P ^ip) l^'ip 



(22) 

where / and / refer to the initial and final value, respectively, 
for the input parameter kp^s2 and output A2p. Using Eq. ( l22l l we 
have calculated the amplification factor as a function of input 
stimulus. The resultant data are plotted in Figure |5]which shows 
a gradual amplification of TF as external stimulus is increased. 
From Figure |5] it is also evident that due to high stimulus the 
gain raises upto ~ 42% compared to its value at low signal. 

3.3. Time dependent dynamics 

To study the time dependent behavior of different quanti- 
ties of the model, the sets of nonlinear ODEs are solved by 
XPP (http://www.math.pitt.edu/~bard/xpp/xpp.html) using the 
parameter set given in Table 2. The parameters listed in Ta- 
ble 2 were guessed to generate the temporal experimental pro- 
file given in Figures (l6] |7] H] and |9^). The consensus set 
of parameters used for numerical integration of the nonlinear 
ODEs are obtained using Parameter Estimation Toolkit (PET) 
(http://mpf.biol.vt.edu/pet/) . 

3.3.1. The bvg operon 

In Figure |6l we compare the nu merical results with exper- 



imental data (IScarlato et al.L 119911) . for time evolution of tran- 
scripts niAs I and mAs2 generated by the two promoters Pas i and 
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Figure 6: Time evolution of transcripts niASl and "Ms 2 due to Pa51 and Pas 2 
activity, respectively. The open symbols aie taken from IScarlato et alj il99lh 
and the solid lines are results of numerical integration. 



Figure 8: Time evolution of four different classes of mRNA due to differe ntial 
gene regulation. The open symbols are taken from IScarlato et all jl99lh and 
the solid lines are due to numerical integration. 



Pa52, respectively. While plotting the data we have scaled each 
transcripts value by the maximum of niAST {- mAsi + mAsi) 
to compare simulated data with experimental results in a rela- 
tive scale of 100. From |6] it is evident that our model captures 
the qualitative aspects of in vivo experimental results. To get an 
idea of how well our model can mimic the real system we define 
the amplification factor, / (= induced/basal) for mRNA expres- 
sion which for experiment and simulation are ^ 2.63 and 



ment. 



2.15, respectively, and are in good agree- 




w 30 h 

CO 

b 20 h 
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3.3.2. The proteins 

As the bvg operon gets switched on at higher temperature 
(37 °C) it starts producing a large pool of response regula- 
tor BvgA which is about ^ 50 fold hig her than that what is 
produced at the non- i nducin g condition (IScarlato et al. . Il990[ 



199 It iPrugnola et al.L 1 19951) . Experimental results show low 



level synthesis of sensor (BvgS) and response regulator (BvgA) 
proteins at low temperature (25 °C), but 6 hours after induc- 
tion their level becomes 56 and 4 fold higher, respectively. 
In the next 18 hours level of BvgA protein goes down to 40 
fold whereas level of BvgS protein finally increases to 17 fold 
dScarlato et al.L Il99ll) . It is important to note that downfall 



of BvgA level after first 6 hours is a signature of arrival of 
B. pertussis colony at the stationary state where nutrition and 
other growth resources may become limited, that effectively 
lead to nonl i near d egradation ( Monodi 1949; Tan et al. . 2009 



Ghosh et al.L 1201 ll) . In the present model we have considered 



only linear degradation kinetics as the relevant dynamics for 
the activation of the downstream genes takes place within first 
8 hours of induction (see discussion in the following). In other 
words, within the first 6-8 hours the pool of BvgA becomes sat- 
urated enough to trigger the signal transduction network. Keep- 
ing this in mind, the present model with linear degradation ki- 
netics can still take care of the relevant dynamics within the ex- 
ponential growth phase and at the transition from exponential 
growth phase to stationary phase. In Figure Q we show quali- 
tative agreement of simulated result with that of experimental 
data for fold increase of total BvgS {Sj) and total BvgA {Aj) 
protein. While caculating St and Aj we have used relations 
(Eq. (14)-Eq. (15)). 



Figure 7: Fold increase of total amount of sensor (Sr) and response regulator 
(A 7-) protein. The dynamic s is shown for first 1 8 hours after induction. The 
open symbols are taken from lScarlato et al1 ( ll99ll) and the solid lines are results 
of numerical integration. 



3.3.3. Gene regulation of phenotypic phases 

After 2 hours of induction total BvgA level increases only 
by ~18 fold, but such low level of response regulator protein 
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is enough to activate class 2 genes that encode the proteins for 
adherence due to high affinity binding site upstream of TSP. As 
a result, within 2 hours of induction the promoter for class 2 
gene gets activated and the corresponding mRNA (111^12) level 
reaches its maximum value (see Figure |8]l. After this 2 hours 
time window, level of BvgA rises and the accumulated amount 
is enough to bind the low affinity binding site of class 1 genes 
and to activate them. Thus, in a period of 2-6 hours of acti- 
vation, class 1 genes gets fully on leading the corresponding 
mRNA (nicn) level to its maximum value (see Figure|8]l. 

Together with class 1 and class 2 genes we have shown the 
time evolution of class 3 and class 4 genes in Figure[8] As men- 
tioned in the introduction, although class 3 gene expression has 
not been observed in B. pertussis under temperature induction, 
we predict that proper modulation of external temperature may 
lead to expression of class 3 gene in B. pertussis as gradual tun- 
ing of external temperature slowly acc umulates the transcrip - 
tion factor in a graded response manner (iPrugnola et ali Il995b . 
Similarly, expression of class 4 gene can be achieved by in- 
corporating a plasmid containing frlAB fused with lacZ from 
B. bronchiseptica into B. pertussis. Similar kind of technique 
has been used t o express ptx promoter from B. pertussis in B. 
bronchiseptica (IWilhams and Gotten. l2007h . 



3.3.4. The mutants 

In the previous discussion we have provided an account of 
how well the developed model can reproduce the qualitative 
features of signal transduction at the molecular level. Next we 
check the vaUdity of the model b y looking in t o soni e novel mu- 
tants reported in the literature (J ones et al.l 120051) . As men- 
tioned earlier necessary condition for the activation of signal 
transduction in B. pertussis is the phosphorylation of the TP, 
BvgA. Thus, any form of hindrance/activation through mu- 
tation at the phosphorylation site of BvgA will get reflected 
in their ability to activate the signaling machinery. At this 
point it is important to mention that in an earlier in vitro set 
up it has bee n observed that BvgA gets phosphate from GST 
tagged BvgS (lUhl and Millen.ll994l) . Keeping this in mind two 
mutants, R152H and T194M along with the wild type BvgA 
have been used to study an in vitro phosphorylation kinetics 
jjones et all l2005h where in a mixture of 0.8 /iM GST-'BvgS 
and 2.1 pM BvgA (wild type or R152H or T194M) 30 ^lU 
[y-^^PJ-ATP was added and the phosphotransfer kinetics was 
monitored at 0.8, 2 and 5 minutes after addition. The relative 
amount of radioactive phosphate incorporated into BvgA was 
then calculated as a fu nction of time usin g phosphoimager (see 
Figures 5B and 5G of IJones et al.l (120051) ) in a relative scale of 
100. The mutant R152H has been reported to behave almost 
like wild type in their ability of getting phosphorylated whereas 
T194M has been shown to be phosphorylated in a drastically 
low amount, ~ 20% of that of wild type BvgA. 

In our model one can control the rate of phosphorylation of 
TF through the Michaelis constant KmA- {kt,b + k,^ai)lktj) for 
the kinase activity of the sensors. For the two mutants R152H 
and T194M we take Kmi to be 3.42 nM and 0.5 nM, respec- 
tively, compared to 9.96 nM for wild type strain. The in vitro 
phosphorylation assay results can be modeled in the present 



study using only Eqs. (I8]|9]), as Eq. ^ mimics phosphoryla- 
tion of GST tagged BvgS using [y-^¥]-ATP and Eq. ^ takes 
care of kinase activity of BvgS towa rds BvgA. To rep roduce 
the in vitro kinetic data reported by Jones et al.l (120051) using 
our model we have used 0.8 pM 5 2 and 2.1 pM. A2 (for wild 
type, R152H and T194M) and varied the Michaelis constant 
Kmi which mimics the role of radioactive ATP incorporation 
(from 5 2 to A2) in the experimental setup. In Figure |9j a) we 
show the in vitro phosphorylation assay results (symbols) along 
with data generated using Eqs.([8]|9]l (lines), which shows a good 
agreement between theory and experiment. 

Nonspecific binding between high affinity binding site and 
BvgA can occur even when the later is unphosphorylatedbut its 
specific affinity for D NA increases in the phosphorylated form 



(iBoucher et al.l 1 19971) . From this information one can expect 



phosphorylated BvgA with R152H substitution will have higher 
binding probability to the primary high affinity BvgA binding 
site, compared to T194M substitution. But surprisingly, almost 
reverse result was observed in ele ctrophoretic mobility shift as- 
says (EMSA) dJones et al.Ll2005h . In the EMSA experiment an 
oligonucleotide (22S YM) with a perfect inverted heptameric re- 
peat at the center was used as it has bee n shown to represent 
a high affinity BvgA binding site earlier ( Boucher and Stibit"3 



119951 



g n aitinity t 
5t IBoucher ( 



et all 119971 1200 lb . Resulting data shows that. 



compared to wild type, BvgA with T194M substitution shows 
~ 40% binding ability compared to wild type. Whereas, BvgA 
with R152H substitution shows almost zero binding affinity to- 
wards the high affinity DNA binding site. The in vitro phospho- 
rylation assay and electrophoretic mobiUty shift assay results 
show reverse effect for R152H and T194M when it comes to 
binding probability to primary high affinity binding site. Thus, 
in vitro results suggests that although replacement of arginine 
(R) by histidine (H) at 152 residue retains the positive charge, 
its interaction with negatively charged double helix backbone 
reduces severely. On the other hand replacement of threonine 
(T) by methionine (M) at 194 position increases the hydropho- 
bicity of the TF so that inspite of being poorly phosphorylated 
their interaction with high affinity binding site increases. 

The above mentioned results with reverse effects thus imme- 
diately raises the question - how these two mutants will behave, 
if expressed in vivo! To answer this we have predicted the tem- 
poral behavior of TF, class 2 mRNA and class 3 mRNA for 
R152H and T194M substitution and compared them with wild 
type profiles (see Figure |3b)-|9td)). While simulating the full 
network (using Eqs. (1-13)) we have decreased all the binding 
and unbinding {kb and fc,,) rates of the model, listed in Table 2, 
by two order of magnitude to get desired phenotypic behavior 
for the mutants R152H and T194M. Our in silico prediction 
suggests that opposing effect of phosphorylation and binding 
affinity causes delay in t he activation of the signal transduction 
machinery as reported in lJones et al.l(l2005l) . 



3.3.5. Sensitivity Analysis 

To check the sensitivity of the parameter values on the net- 
work dynam ics listed in Table 2 we hav e adopted the procedure 
described bv lBarkai and Leiblerl(ll997l) . In this procedure all or 
a subset of the rate constants were subjected to a random pertur- 



8 



I I I I I 

(a) WT — 

R152H — 
T194M ■■■■ 



_i I I I i_ 



WT - 
R152H ■ 
T194M ■ 



time (min) 
1 1 1 r- 



1 2 3 4 5 6 
time (hr) 






100 


. ( 






















80 




/ \ WT 








' \ R152H 




60 




\ T194M 


E 










40 






> 








IT 


20 
















12 3 4 5 6 
lime (hr) 



CL 
CM 

< 



247 p 

246.5 - 

246 - 

245.5 - 

245 - 

244.5 - 

244 - 
243.5 



2 4 6 8 10 12 14 16 18 20 
perturbation (%) 



Figure 9: Effect of mutation on pliosphorylation of BvgA and its effect on tlie 
activation of genes with liigli affinity bi nding site, (a) Compaiison of in vitro 
pfiosplioryiation data (symbois are from jjones et al.l ll2005)) and theoreticai re- 
suits (soiid, dashed and dotted iine) for WT, R152H and T194M. (b) Reiative 
amount of A2p due to WT and the two mutants R152H and Ti94M. (c) and 
(d) Reiative amount of transciipts generated from genes with high affinity bind- 
ing site of TF. In all the four panels all the data are scaled with respect to the 
maximum value, set to 100. 



bation where the perturbation was drawn from a random gaus- 
sian distribution whose mean is the unperturbed value of each 
rate constant and variance is certain percentage (up to maxi- 
mum of 10%) of the rate constant. This leads to a set of unper- 
turbed (reference) rate constants, A:" (listed in Table 2) and an 
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Figure 10: Plot of steady state [A2p] as a function of total parameter variation 
\og(k). The horizontal dashed line is the steady state A2p value (~ 244 nM) for 
the reference system whereas each dot represents the same for modified set. (a) 
Only binding constants are modified; (b) only synthesis and degradation rates 
are modified; (c) only kinase and phosphatase rates are modified and (d) all the 
rates are modified. 



Figure 1 1 : Evolution of mean [A2p] level as function of perturbation (%) (dot- 
ted line). The horizontal solid line is the steady state A2P value (~ 244 nM) for 
the unperturbed system. 



ensemble of perturbed (modified) set of rate constants, A:,. We 
then calculated the level of A2P at steady state using both the 
reference set and the modified set using the full network. Varia- 
tion of the reference steady state [A2p] value for the modified set 
of rate constants can be characterized by total parameter varia- 
tion k , defined as log(fc) - Yjf=i |log(A;,7^")| (Barkai and Leibler, 
1997h . The resultant simulation results are showed in Figure [TOl 
where the dotted line is the steady state [A2p] for the reference 
state and each dot represents the same using the perturbed set. 
We have tested the sensitivity of the parameters using four dif- 
ferent sets. In set 1 we have only perturbed the binding con- 
stants (Figure [TOk); in set 2 the rate constants controlling the 
synthesis and the degradation of the systems components have 
been modified (FigurefTOb): in set 3 the rate constants for kinase 
and phosphatase activities have been changed (FigurefTOb) and 
finally in set 4 all the rate constants given in Table 2 have been 
modified (FigurefTOb). From the sensitivity test it is evident that 
the rate constants responsible for kinase and phosphatase activ- 
ities are most sensitive to random perturbation (see the range of 
ordinate in FigurefTOb) compared to the other parameters of the 
model. 



To check further whether the parameter set given in Table 2 
can withstand larger perturbation and gives physically realiz- 
able quantity we have perturbed all the parameters up to 20%. 
Using the perturbation method mentioned before we have gen- 
erated 1000 data set for each level of perturbation and calcu- 
lated the mean [A2p] level to compare with the base value. In 
Figure ( fTTT i we show the evolution of mean [A2p] level as a 
function of perturbation, which shows that even with perturba- 
tion as large as 20% the mean protein level (dotted line) lies 
within the 1% of the unperturbed value (the solid horizontal 
line). Beyond 20% perturbation the solution diverges by sev- 
eral orders of magnitude and becomes unrealistic. 
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4. Conclusion 



5. Appendix 



The kinetic model developed in the present study takes care 
of signal transduction and phenotypic gene regulation in B. per- 
tussis under the influence of temperature elevation. The pro- 
posed model includes all possible elementary kinetics of the 
biochemical interactions between several system components. 
To understand the molecular switch operative in the BvgAS 
TCS, a quasi steady state analysis has been performed which 
reveals temperature induced sharp molecular switch that re- 
sponds to the external stimulus, a reminiscent of amplified sen - 



sitivity (iGoldbeter and KoshlandLllQSlUKoshland et alill982h 



Development of the sharp switch has been shown to be the con- 
sequence of positive feedback motif present within the bvgAS 
operon that becomes operative when the temperature of the sur- 
rounding is increased (Figure |5). Outcome of the sharp switch 
gets reflected in the accumulation of the TF in large amount 
within few hours of induction (Figure |7]i, which then controls 
the expression of several downstream genes including the genes 
for adherence and toxin. All these features have been observed 
via numerical integration of the coupled nonlinear ODEs for 
large set of parameters. The resultant numerical results essen- 
tially capture the qualitative features of the network dynamics 
performed in vivo. On the basis of our developed model we 
then looked into the behavior of two novel mutants impaired in 
their ability to phosphorylate the transcription factor and made 
testable predictions for temperature induced class 3 gene ex- 
pression. 

We hope that our in silico study will inspire more experi- 
ments in the coming days to address subtle issues of the network 
that are yet to be explored. One of such issues is the character- 
ization of exact target and functioning of the antisense RNA, 
transcribed from Pasa promoter of bvgAS operon. In this con- 
nectio n it is important to mention the recent work by iHot et al 



(1201 lb . where sRNA of B. pertussis have been identified for the 
first time. Out of the pool of sRNA identified one ibprJI) has 
been mentioned to be controlled by BvgAS TCS but exact tar- 
get and mode of its functionality is yet to be discovered. In 
addition, one may find it interesting to make quantitative mea- 
surement of different proteins generated due to activity of the 
four classes of downstream genes. Information from these new 
experimental findings will certainly help one to build a more 
detailed model in future. 
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5.7. Promoter kinetics of downstream genes 

The active and inactive forms of four different promoters 
of the downstream genes are modeled as Pcij,a and Pdj.i ij - 
1,2,3,4), respectively. Considering co-operativity in binding 
of TF (A2/)) to the high/low affinity binding site of these pro- 
moters we model binding kinetics as follows. 
For class 1 gene: 



PcllJ + Mp ^ Pd\,i\, 

Km 

h.ii 

Pclljl + Mp ^ Pclljl, 



Pclljl + Mp ^ Pcn,a- 



For class 2 gene: 

Pcll,i+A2P ''^ Pclljl, 

ki,.i2 

Pclljl + Mp ^ Pclljl, 



Pclljl + AiP ^ Pcll,a- 



For class 3 gene: 



PcBJ+Aip ^ PcBjl, 

h.i2 

PcBJl + Aip ^ Pcn,a, 

k„.32 

kb.3i 

PcB,a + A-IP ^ PcBJl, 



For class 4 gene: 



Pcl4,a +Aip^ PclAJ, 



(23a) 
(23b) 
(23c) 

(24a) 
(24b) 
(24c) 

(25a) 
(25b) 
(25c) 

(26) 



In the above set of equations Pdjjk {k - 1,2,...) represents the 
inactive intermediate states of the different classes of promot- 
ers. 

5.2. Promoter kinetics ofbvg locus 

As mentioned in the main text we consider a single copy of 
file bvg gene, with relations [Pasic] + [Pasic] + [Pash] = 1 
and [Pasu] + [^'asio] - 1 for the two promoters controlling 
functionality of the bvg operon. Dynamical equations for the 
promoter kinetics thus can be written as 



djPASlc] 

dt 

d[PASla\ 

dt 



d[PASla\ 

dt 



-hl[PASlc][Alp]+k,i[PASla], (27) 
hl[PASlc][Alp] + kal[PASli\ 

-(kui + ka)[PAsia], (28) 

kbl[PASu][Alp]-k,APASla], (29) 
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which at the steady state give, 

[Pas la] = Fi(A2p), [Pasu] = 1 - i^i(A2/>), 

[PaS2c] = F2(A2P), [PASla] = i^2(A2p)[A2p]2, (30) 

where 

1 



^i(A2f) = . ,F2(A2p) 

1 + A2P 1 



l+(l+k)[A2p]2 



with [A2p]i = [A2P]/^:; for Ki = k,i/kbi (i = 1,2) and A: = 
^;2/^fl2- 

5.3. niRNA kinetics 

Time evolution of the transcripts generated from the two pro- 
moters Pasi and Pas I are given by the following set of equa- 
tions, respectively, 

diniAsi] 



dt 

d[mAs\] 
dt 



- k,p^20[PAS2c] + k,p^2l[PAS2a\ 
-kd,m[mAS2], 

- k,p^ 1 1 [Pas \a\- kd.m [mAs i ] ■ 



(31) 
(32) 



Along with the steady state expression for the two promot- 
ers given in Eq. ( l30t . one arrives at the following expres- 
sions for concentrations of the two transcripts at steady state 
id[mAS2]/dt - d[mAsi]/dt - 0) 



= -jf F2iA2p) + -f i^2(A2p)[A2p]2, 

ktp 11 

[mAs\] = -j-^F\{A2p). 



(33) 



5.4. Phosphotransfer kinetics 

From the kinetics of phosphate donation from sensor to re- 
sponse regulator (kinase), and phosphate withdrawal from re- 
sponse regulator by sensor (phosphatase), we have two dynam- 
ical equations for the two intermediates S2P ■ A2 and 5 2 ■ A2P 



d{S2p-A2\ 
dt 

d[S2-A2p] 
dt 



= k,j{S2p]{A2\ 

-{k,^b+k,^a2)[S2P ■ M], 

= kpj[S2][A2p] 

-(kp^b + kp,a2)[S 2 ■A2p]. 



(34) 



(35) 



While writing the above two equations we have neglected 
degradation of the intermediates as they show transient dy- 
namics, i.e., time evolution of both the Michalies complexes 
takes place on a faster time scale compared to the other system 
components. Now by imposing quasi-steady state conditions 
d[S2p-A2]/dt = and d[S2-A2p]/dt = on Eqs. ( I34ll35l l we 
have 

, , [^2f][A2] , . [^2][A2P] 

[S2P-A2\- — ,li2-A2pJ = — , (36) 



K, 



Mp 



where Kmi = (ktb + k,^a2)lk,j and Kmp = {kp^b + kp,a2)/kpj 
are the Michaelis constants for the kinase and the phosphatase 
activity of the sensors, respectively. 



5.5. Protein kinetics 

The dynamical equations representing kinetics of the differ- 
ent forms of sensor and response regulator proteins can be rep- 
resented by the following sets of ordinary differential equations 



d[S2] 
dt 



d[S2p] 
dt 

d[A2] 
dt 



d[A2p] 
dt 

where 

a I = ^.5.5,2 



aiF2(A2p) + Q'2f2(A2/')[A2p]2 + aiFi(A2p) 



kt. 



,£l2 , 



-kp,s2{s)[S2]+kdp,s2[S2p] + -i^[S2p][A2\ 



-kd,p[S2], 



(37) 



kL6 



kp,s2i^)[S2] -kdp,s2[S2p] - ■7^[S2P][A2] 

J^Mt 

-kd.p[S2p], (38) 

jei^2(A2p) +^2i^2(A2F)[A2p]2 +kFl(A2p) 

-tF[52p][A2] + ^[S2\[A2P] 
J^Mt J^Mp 



-kd,p[A2], 

-^{S2P]{A2\ - -^{S2\{A2P] 

-kd,p[A2p], 



Mp 



(39) 



(40) 



klp,2{) _ , ^/;),21 . , k,p^i\ 

— ,a2 - kss,2^ ,^3 - «i.!,lT ' 

'^d,m '^d,m t^d,!)! 

ktp,ll 



3 ktp,20 ~ k,p^2\ ~ 

Pi - K'sa,2~r 'P2 - K-sa,2— ,P3 - K.s£i,l , 

t^d^m '^d,m '^d,m 

While writing Eqs. (l37H40b we have used the steady state ex- 
pressions for [mAs\], ['«AS2], [S2P ■ A2] and [^2 • A2p] given by 
Eq. ^ and Eq. 



5.6. Linear stability analysis 

Using the relations [^2] + [S2p] ~ [Sp], [A2] + [A2p] ~ [Ap] 
and Eqs. (I37H40I I one can write 



djSp] 

dt 
djAp] 

dt 



f(A2p)-kd,p[ST]^0, 
g{A2p) - kd,p[Ap] = 0, 



(41) 
(42) 



where 



/(A2p) = aiF2(A2p) + a2^'2(A2p)[A2p]2 
+a-iFi{A2p), 

8(A2P) = j§l/^2(A2p)+y62F2(A2p)[A2p]2 
+j63^'l(A2p). 



(43) 



(44) 



At 25 °C when the switch is off, A2P - and hence f{A2p) = 
ai + a3, and ^(A2p) = $1 +^3,- Thus we have at steady state the 
fixed point 

[5r].„ = ai, [At]ss =y8i 
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where a, - &i/k4^p and yS; - ^i/kd,p for 1 - 3. On the other 
hand, at 37 °C for A^p, the steady state value of A2p, we have 

To understand the nature of the fixed point ([StTss> IAtYss) we 
construct the stability matrix evaluated at steady state, 

-kd,p \ 
-kd,p j' 

for which trace, t = -2kd^p{< 0) and determinant, A = ^p(> 
0), a characteristics of stable fixed point. 
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Table 1 : List of symbols (with initial values) used in the model 



Symbol 


Initial Value 


Description 


Pas la 


OnM 


Active state of promoter Pas 1 


Pasu 


0.98 nM 


Inactive state of promoter Pas 1 


Pas 2c 


0.98 nM 


Constitutive state of promoter Pas2 


Pas 2a 


nM 


Active state of promoter Pas2 


PAS2i 


OnM 


Inactive state of promoter Pas2 


niAsi 


OnM 


Transcripts generated from Pas la 


mAS2 


1.11 nM 


Transcripts generated from Pas 2a 


S2 


10.74 nM 


Dimer of sensor BvgS 


Ai 


11.23 nM 


Dimer of response regulator BvgA 


S2P 


nM 


Phosphorylated dimer of sensor BvgS 


Alp 


OnM 


Phosphorylated dimer of response 






regulator BvgA (transcription factor) 


S2P-A2 


nM 


Michaelis complex formed by S2P and A2 


S2 ■ A2P 


OnM 


Michaelis complex formed by 5 2 and Aip 


Pen,, 


OnM 


Active promoter of class 1 gene 


PcllJ 


0.98 nM 


Inactive promoter of class 1 gene 


Pcl2,a 


OnM 


Active promoter of class 2 gene 


Pcl2,i 


0.98 nM 


Inactive promoter of class 2 gene 


PcB,a 


OnM 


Active promoter of class 3 gene 


PcB,i 


0.98 nM 


Inactive promoter of class 3 gene 


Pcl4,a 


0.98 nM 


Active promoter of class 4 gene 


Pcl4,i 


OnM 


Inactive promoter of class 4 gene 


men 


OnM 


Transcripts generated from Pcii,a 


mci2 


OnM 


Transcripts generated from Pci2,a 


mcB 


OnM 


Transcripts generated from PcB,a 


rriciA 


2.93 nM 


Transcripts generated from Pc/4,a 
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Table 2: List of kinetic parameters (with values) used in the model 



Parameter 


Value 






Description 


kb\ 


1.024 X 10" 


"71 TZ — 

* nM" 


n 1 — 

' s"' 


Association rate of A2p and Pas u 


ku\ 


1.167 X 10" 


"3 S^^ 




Dissociation rate of A2p from Pasu 


hi 


1.36 X 10"^ 


nM-i 


s-1 


Association rate of A2p and Pas2c 




2.5 X f 






Dissociation rate of A2p from PAS2a 


ka 


1.667 X 10" 


3 s-i 




Inactivation of P2a to f as2i 




2.0 X 10""* f 






Activation of P2a from PAS2i 


ktp,20 


1.9 X 10"^ f 


1 

1 — 1 




Transcription rate from Pas2c promoter 


ktp,2\ 


9.386 X 10" 


^ nM" 


' s"' 


Transcription rate from PAS2a promoter 


ktp,n 


4.083 X 10" 


S-' 




Transcription rate from Pasxo promoter 


kd,m 


1.667 X 10" 


-3 S-' 




Degradation rate of mRNA 


kss,l 


6.667 X 10" 


S-' 




Synthesis rate of 5 2 from otasi 


kss,2 


1.667X 10" 


-3 s-^ 




Synthesis rate of 5 2 from mAS2 


ksa, 1 


4.167 X IQ- 


■2 S-' 




Synthesis rate of A2 from mas 1 


ksa,2 


1.667 X 10" 


"3 s"' 




Synthesis rate of A2 from mAS2 


kp,s2{s) 


8.333 X 10" 


-3 s"^ 




Phosphorylation rate of 5 2 at 37 °C 


kdp,s2 


3.333 X 10- 


-3 S-' 




Dephosphorylation rate of S2P 


ktj 


8.532 X 10" 


3 nM" 


' S-' 


Association rate of S2P and A2 


kt,b 


1.667 X IQ- 


-3 s-i 




Dissociation rate of 52/' • A2 


kt,a2 


8.333 X IQ- 


■2 s-i 




Phosphate transfer rate from S2P to A2 


kpj 


3.413 X 10" 


5 nM- 


1 S-' 


Association rate of 5 2 and A2P 


kp,b 


1.333 X IQ- 


-3 s-i 




Dissociation rate of ^2 • A2P 


kp,a2 


5.0 X 10-2 s 


,-1 




Phosphate removal rate from A2P by 5 2 


kd,p 


1.667 X 10" 


"4 s-i 




Degradation rate of protein 


kb,n 


6.826 X 10" 


-^nM- 


^ s-1 


Association rate of A2P and Pdu 


ku,n 


1.667X 10" 


-6 ^-1 




Dissociation rate of A2P from Pcii,a 


kb,\2 


1.024 X 10" 


nM- 


1 s-1 


Association rate of A2P and Pcn,n 


ku,\2 


1.667 X 10- 


■6 ^-1 




Dissociation rate of A2P from Pcn,i2 


kb,n 


1.36 X 10"^ 


nM-i 


s-1 


Association rate of A2P and Pd\,i2 


ku,n 


1.667 X 10" 


-6 g-1 




Dissociation rate of A2P from Pci\,a 


kh,2\ 


5.119X 10" 


"'nM- 


s"^ 


Association rate of A2P and Pci2,i 


ku,21 


1.667 X 10- 


-4^-1 




Dissociation rate of A2P from Pci2,ii 


kb,22 


1.36 X 10^ 


nM-i 


s-1 


Association rate of A2P and Pci2,n 


ku,22 


1.667X 10- 


-4 ^-1 




Dissociation rate of A2P from Pci2,i2 


kb,23 


1.706X 10" 


■^nM- 


1 s-1 


Association rate of A2P and Pci2,i2 


ku,23 


1.667 X 10" 


S-' 




Dissociation rate of A2P from Pci2,a 


kh,3i 


8.533 X 10" 


nM" 


1 1 

^ s"^ 


Association rate of A2P and PcB,i 


ku,3l 


1.667 X 10- 


A 1 




Dissociation rate of A2P from Pci3,i\ 


kb,32 


1.365 X 10" 


^ nM" 


1 1 

^ s"^ 


Association rate of A2P and Pci3,n 


ku,32 


1.667 X 10" 


^ s ^ 




Dissociation rate of A2P from Pci3,a 


kb,33 


1.706 X 10- 


■^nM- 


1 s-1 


Association rate of A2P and PcB,a 


ku,33 


2.0 X 10-"* f 


,-1 




Dissociation rate of A2P from Pci3,a 


kbM 


1.706 X 10" 


"♦nM- 


1 s-1 


Association rate of A2P and PcH,a 


kuM 


1.667X 10- 


-4^-1 




Dissociation rate of A2P from f cM,/ 


ktp,cn 


5.167 X 10" 


-3 s-i 




Transcription rate from Pci\,a promoter 


ktp,cl2 


5.083 X 10" 


-3 s-i 




Transcription rate from Pci2,a promoter 


ktp,cl3 


6.16 X 10-3 


s-1 




Transcription rate from Pci3,a promoter 


ktp,clA 


5.0 X 10-3 s 


,-1 




Transcription rate from PcIA^ promoter 
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